library(tidyverse)  # ggplot(), %>%, mutate(), and friends 
library(broom)
library(MatchIt)  # Match things
library(Rcpp)
library(MASS)
library(modelsummary)
library(IRdisplay)
library(haven)
library(marginaleffects)
library(dplyr)
library(jtools)
library(ggeffects)
library(jtools)
library(lme4)

# set directory


# IMPORT datasets. equal and discordant
equal <- read_csv('C:/./dataset_map_equal_ready_with_ses.csv')
disc <- read_csv('C:/./dataset_map_discordant_ready_with_ses.csv')


################################
################################
################################
# E Q U A L  M A P  A N D  W P #
################################
################################
################################
library(janitor)

#clean labels: can be done by simply
equal <- clean_names(equal)




# check whether means are statistically significant via t-test
with(equal, t.test(solved_yes ~ vic_race_black))


# exposure as character
equal$vic_race_black <- as.character(equal$vic_race_black)



equal2 <- equal %>% drop_na()
#################################################################
####################### MAIN MODELS #############################
#################################################################




############## MAP (agency info + states) (adjusted, unmatched)
m.out3wp <- matchit(vic_race_black ~ vic_count + 
                    off_count + 
                    victim_sex+
                    age_cat+
                    decade+
                    agentype+
                    monthly_agency_overlap+
                    state,
                  data=equal,
                  method = "exact",
                  estimand="ATT")

bal_tab_3wp<-bal.tab(m.out3wp, un=TRUE)


matched_dataset3wp <- match.data(m.out3wp)



# Adjusted (with ses)
model_matched3wps <- glm(solved_yes ~ vic_race_black + 
                        vic_count + off_count + victim_sex+
                        age_cat + 
                        decade+
                        weapon+
                        circumstance+
                        #agentype removed because has only one level
                        monthly_agency_overlap+
                        state+
                        acs_income_higher100+
                        acs_demographics_race_and_ethnicity_not_hispanic_or_latino_black_or_african_american_alone_percentage,
                      family=quasibinomial(), data = matched_dataset3wp, weights=weights)

tidy(model_matched3wp)



summ(model_matched3wps, exp = TRUE, robust=TRUE, digits=3)
ame3wps <- comparisons(model_matched3wps, variables = "vic_race_black", vcov="HC3",
            wts="weights",
            newdata = subset(matched_dataset3wp, vic_race_black == 1)) |>
  summary()

ame3wps$dataset <- 'MAP-WP (Equal Outcomes) w/Census'
ame3wps$model <- 'Adjusted/Matched'

# Adjusted (no ses)
model_matched3wp <- glm(solved_yes ~ vic_race_black + 
                           vic_count + off_count + victim_sex+
                           age_cat + 
                           decade+
                           weapon+
                           circumstance+
                           #agentype removed because has only one level
                           monthly_agency_overlap+
                           state,
                        family=quasibinomial(), data = matched_dataset3wp, weights=weights)

summ(model_matched3wp, exp = TRUE, robust=TRUE, digits=3)
ame3wp <- comparisons(model_matched3wp, variables = "vic_race_black", vcov="HC3",
                       wts="weights",
                       newdata = subset(matched_dataset3wp, vic_race_black == 1)) |>
  summary()

ame3wp$dataset <- 'MAP-WP (Equal Outcomes) no Census'
ame3wp$model <- 'Adjusted/Matched'


# Unmatched with ses
model_unmatched3wps <- glm(solved_yes ~ vic_race_black + 
                          vic_count + off_count + victim_sex+
                          age_cat + 
                          decade+
                          weapon+
                          circumstance+
                          monthly_agency_overlap+
                          state +
                          acs_income_higher100+
                          acs_demographics_race_and_ethnicity_not_hispanic_or_latino_black_or_african_american_alone_percentage, 
                        family=binomial, data = equal)


summ(model_unmatched3wps, exp=TRUE, robust=TRUE, digits=3)
ame3wps_u <- comparisons(model_unmatched3wps, variables = "vic_race_black", newdata = subset(equal, vic_race_black == 1)) |>
  summary()

ame3wps_u$dataset <- 'MAP-WP (Equal Outcomes) w/Census'
ame3wps_u$model <- 'Adjusted/Unmatched'


# Unmatched no ses
model_unmatched3wp <- glm(solved_yes ~ vic_race_black + 
                            vic_count + off_count + victim_sex+
                            age_cat + 
                            decade+
                            weapon+
                            circumstance+
                            monthly_agency_overlap+
                            state,
                           family=binomial, data = equal)


summ(model_unmatched3wp, exp=TRUE, robust=TRUE, digits=3)
ame3wp_u <- comparisons(model_unmatched3wp, variables = "vic_race_black", newdata = subset(equal, vic_race_black == 1)) |>
  summary()

ame3wp_u$dataset <- 'MAP-WP (Equal Outcomes) no Census'
ame3wp_u$model <- 'Adjusted/Unmatched'

################################
################################
################################
###### D I S C O R D A N T #####
################################
################################
################################

#clean labels: can be done by simply
disc <- clean_names(disc)




# check whether means are statistically significant via t-test
with(disc, t.test(solved_yes ~ vic_race_black))


# exposure as character
disc$vic_race_black <- as.character(disc$vic_race_black)



disc2 <- disc %>% drop_na()

#################################################################
####################### MAIN MODELS #############################
#################################################################



############## MAP (agency info + states) / Unadjusted, adjusted, unmatched
m.out3dwp <- matchit(vic_race_black ~ vic_count + off_count + victim_sex+
                    age_cat+
                    decade+
                    agentype+
                    monthly_agency_overlap+
                    state,                        
                  data=disc,
                  method = "exact",
                  estimand="ATT")


bal.tab(m.out3dwp)

matched_dataset3dwp <- match.data(m.out3dwp)




# Adjusted with ses
model_matched3dwps <- glm(solved_yes ~ vic_race_black + 
                        vic_count + off_count + victim_sex+
                        age_cat + 
                        decade+
                        weapon+
                        circumstance+
                        #agentype	 (only one level, removed)
                        monthly_agency_overlap+
                        state+
                        acs_income_higher100+
                        acs_demographics_race_and_ethnicity_not_hispanic_or_latino_black_or_african_american_alone_percentage, 
                      family=quasibinomial(), data = matched_dataset3dwp, weights=weights)

tidy(model_matched3dwps)

summ(model_matched3dwps, exp = TRUE, robust=TRUE, digits=3)
ame3_dwps <- comparisons(model_matched3dwps, variables = "vic_race_black",
                          wts="weights", newdata = subset(matched_dataset3dwp, vic_race_black == 1)) |>
  summary()

ame3_dwps$dataset <- 'MAP-WP (Discordant Outcomes) w/Census'
ame3_dwps$model <- 'Adjusted/Matched'

# adjusted no ses
model_matched3dwp <- glm(solved_yes ~ vic_race_black + 
                           vic_count + off_count + victim_sex+
                           age_cat + 
                           decade+
                           weapon+
                           circumstance+
                           #agentype	 (only one level, removed)
                           monthly_agency_overlap+
                           state,
                         family=quasibinomial(), data = matched_dataset3dwp, weights=weights)

tidy(model_matched3dwp)

summ(model_matched3dwp, exp = TRUE, robust=TRUE, digits=3)
ame3_dwp <- comparisons(model_matched3dwp, variables = "vic_race_black",
                        wts="weights", newdata = subset(matched_dataset3dwp, vic_race_black == 1)) |>
  summary()

ame3_dwp$dataset <- 'MAP-WP (Discordant Outcomes) no Census'
ame3_dwp$model <- 'Adjusted/Matched'

# Unmatched with ses

model_unmatched3dwps <- glm(solved_yes ~ vic_race_black + 
                           vic_count + off_count + victim_sex+
                           age_cat + 
                           decade+
                           weapon+
                           circumstance+
                           agentype+	
                           monthly_agency_overlap+
                           state+ acs_income_higher100+
                          acs_demographics_race_and_ethnicity_not_hispanic_or_latino_black_or_african_american_alone_percentage, 
                         family=binomial, data = disc)

summ(model_unmatched3dwps, exp = TRUE, robust=TRUE, digits=3)
ame3_dwps_un  <- comparisons(model_unmatched3dwps, variables = "vic_race_black", newdata = subset(disc, vic_race_black == 1)) |>
  summary()

ame3_dwps_un$dataset <- 'MAP-WP (Discordant Outcomes) w/Census'
ame3_dwps_un$model <- 'Adjusted/Unmatched'


# Unmatched no ses
model_unmatched3dwp <- glm(solved_yes ~ vic_race_black + 
                             vic_count + off_count + victim_sex+
                             age_cat + 
                             decade+
                             weapon+
                             circumstance+
                             agentype+	
                             monthly_agency_overlap,
                           family=binomial, data = disc)

summ(model_unmatched3dwp, exp = TRUE, robust=TRUE, digits=3)
ame3_dwp_un  <- comparisons(model_unmatched3dwp, variables = "vic_race_black", newdata = subset(disc, vic_race_black == 1)) |>
  summary()

ame3_dwp_un$dataset <- 'MAP-WP (Discordant Outcomes) no Census'
ame3_dwp_un$model <- 'Adjusted/Unmatched'

# plotting (Figure on national AMEs)

ame_list <- do.call("rbind", list(ame3wps, ame3wp, ame3wps_u, ame3wp_u,
                                  ame3_dwps, ame3_dwp, ame3_dwps_un, ame3_dwp_un ))
ame_list <- ame_list[order(ame_list$dataset),]
ame_list$dataset<- fct_inorder(ame_list$dataset)
ame_list$model = factor(ame_list$model, levels=c("Adjusted/Matched", "Adjusted/Unmatched"))
# plot
dodge <- position_dodge(width=0.5)
ame_plot <- ggplot(ame_list, x=factor(ame_list$dataset, levels=unique(ame_list$dataset)), y=estimate, color=model, aes(ymin=-0.06,ymax=0.06))+
  geom_point(aes(x=factor(ame_list$dataset, levels=unique(ame_list$dataset)), y=estimate,shape=model, color=model), size=2.5, position=dodge)+
  scale_y_continuous(breaks =c("-.06"=-0.06, "-.04"=-0.04,"-.02"=-0.02, "0"=0, ".02"=0.02, ".04"=0.04, ".06"=0.06))+
  geom_errorbar(aes(x=dataset, ymin=conf.low, ymax=conf.high, group=model, color=model), width=0.1,position=dodge, name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched"))+
  scale_color_manual(values=c("Adjusted/Matched"="red","Adjusted/Unmatched"="blue"), name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched"))+
  scale_shape_manual(values=c("Adjusted/Matched"=19,"Adjusted/Unmatched"= 15), name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched" ))+
  labs(x= "Dataset", y = "Average Marginal Effect (95% C.I.)")+
  theme_bw()+
  geom_hline(aes(yintercept=0), colour="#990000", linetype="dashed")+
  coord_flip()+theme(axis.text.y = element_text(size=12)) + 
  theme(axis.text.x = element_text(size=12), axis.title=element_text(size=16))+
  #theme(text = element_text(family = "Times New Roman"))+
  theme(legend.title=element_text(size=14), legend.key.size=unit(0.5, 'cm'), legend.text = element_text(size=14))+theme(legend.position="bottom")+
  guides(colour=guide_legend(nrow=3,byrow=TRUE))




theme(plot.margin = margin(0.2,0.5,0.2,0.5, "cm"))+theme(legend.position="bottom")

